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Abstract 

We study the nonequilibrium steady state (NESS) in a quantum system in contact with two 
heat baths at different temperatures. We use a time-independent perturbative expansion with re- 
spect to the coupling with the two heat baths to obtain the density matrix for the NESS. In par- 
ticular, we show an explicit representation of the density matrix for the reflection symmetric and 
weakly nonequilibrium case. We also calculate the expectation value of the energy current and 
show that the Kubo formula holds in this case. 



1 Introduction 
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Construction of nonequilibrium statistical mechanics is a challenging open problem in physics. Since 
nonequilibrium phenomena are so diverse, probably it is impossible to make a theory which explains 
all nonequilibrium phenomena. Then a natural first step would be statistical mechanics for nonequi- 
librium steady states (NESSs). The most ambitious goal in this direction is to find a simple theoretical 

O 

o 



expression for the density matrix of the NESSs. 



Although there are many theoretical frameworks to treat NESSs, it is quite difficult to write down the 
density matrix explicitly. For example, in the linear response theory [1] the density matrix for the 
NESS is obtained in the long-time limit of the dynamics under an external field 

oo : t 

Pness = p cq +lim [ dt'e-^^lH^p^-^, (1) 

O i 

and the system has to be infinitely large. (Otherwise we obtain another equilibrium state.) Although 
this equation is useful to calculate some nonequilibrium properties like transport coefficients, Pness 
itself is very hard to calculate in this formalism. Most of formalisms to treat NESS contain this kind 
of long time evolution, which makes it difficult to calculate Pness- 



In this paper, we consider a system with two heat baths. We consider the stationary solution of a 
quantum master equation, and calculate it explicitly using a perturbative expansion with respect to the 
coupling parameter between the system and the heat baths. In particular, in the reflection symmetric 
case we show an explicit form of density matrix for the NESS in the weakly nonequilibrium regime. 



2 Equation of Motion 

We start with the equation of motion for the total system: 

3lAot(t) = T^[H to t, Ptot(t)], (2) 

where 

H tot = Hs + Hb + uH bs . (3) 

Here Hs, Hb and H bs are the Hamiltonians of the system, the heat baths and the interactions, re- 
spectively. We use u as the perturbation parameter. In this paper, we consider a system with two heat 
baths: 

H B = H L + H R , (4) 

Hbs = H LS + H RS - (5) 

Here, Hl and H R are the Hamiltonians for the left and right reservoirs, respectively. We assume 
that the heat bath a {a = L, R) is in equilibrium with the inverse temperature /3 a (See Fig. [[])■ The 
interaction Hamiltonians H LS and Hrs can be written in the form 

H L s = Y.Xft L , H RS = J2X?Y« (6) 

3 3 

where X? acts on the system, and Y L (Y R ) acts on the left (right) heat bath. In the following we 

assume X" and Y-* are Hermitian for simplicity. However, our main results in this paper hold without 
this assumption. 
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Figure 1: The system is in contact with two heat baths at different inverse temperatures L and (3r. 

We expand the density matrix up to 0(u 2 ), trace out the heat bath variables, and apply the Markov 
approximation. Then we obtain the equation of motion for the system [2] 

j t p(t) = ^[H' s ,p(t)}+u 2 J2 ra P(t)- ^ 

a=L,R 



Here, 

H' s = H s + uJ2 E*iW> ( g ) 

a=L,R j 

is the system Hamiltonian with the averaged interaction terms, where (A a ) represents the average of 
A a with respect to the heat bath a. Hereafter H' s is denoted as H s for simplicity. T a is the heat bath 
superoperator, whose explicit form is 



1 f°° < 

r ^(*) = -«E dt'{xfxr{-trp{m%t')-x-p{t)x-{-t')^{-t') 

3,1 J ° 

+p{t)X?{-t')Xf<S>«{-t') - Xn-t'Yp^X^it')) (9) 

Here, X(t) = e~ lHst l h Xe %Ust l h represents the operator in the interaction picture. 

*$(«) = (AY?(t)AY?) (10) 

is a correlation function in the heat bath a, where 

AY* = Yf-(Yf). (11) 

Note that we used 

(AY?(t)AYf) = MA17(£)A17*> (12) 

to derive Eq. ©. Since we assumed Y^'s are Hermitian, 

$ j7 (t)* = *«H). (13) 

The heat bath superoperator © can be rewritten as 

r° = r° + r^, (14) 

where 

r ^ = -^Efft^+ft^)' (15) 

3,1 

T *P = -^Ti^W^-iX^W^). (16) 

3,1 

The operators jfe and W^ are defined as 

(£ P |i$|£,> = (EblXriEjQuM, (17) 

(£ p |W^|£ g ) = (ErlXftEjVfliun), (18) 

where |.E7j) is an energy eigenvector of the system with eigenenergy E*, u vq = (-E p — E q )/K and 



$?,(w) = / rfte- iw4 $"(t) (19) 



oo 



«.(„) = P f ^M^. (20) 



Here V denotes the Cauchy principal value. Note that 

$»* = S«( w ), (21) 

«?,(«)• = *§(«)■ (22) 

The correlation functions satisfy the Kubo-Martin-Schwinger (KMS) condition 

$?,(-") = e^SgM, (23) 

which is equivalent to the following operator identity: 

Rf t = e^ s R%e-^ s . (24) 

Using this identity it is easy to show that 

T a e -/3 a Hs = Q^ 

which guarantees the existence of the equilibrium solution for Eq. © when (3 L = f3 R . 

3 Perturbative expansion 

We put p = const, in the equation of motion ©. Then we have the equation for the steady state 

C p + vfrp = 0, (25) 

where 

C 0P = ^{Hs,p}, (26) 

C 1P = T L p + T R p, (27) 

and v = u 2 . We expand p with respect to v. 

P = Po + vp! + v 2 p 2 + ... (28) 

Then we obtain a series of equations 

£ p = 0, (29) 

C pi + C 1Po = 0, (30) 

C0P2 + C1P1 = 0, (31) 

: (32) 



3.1 Separation of diagonal and off-diagonal parts 

In a normal perturbation theory, we can determine pi step by step starting from the Oth order solution 
p . In this case, however, the Oth order equation (1291 is degenerate, and any diagonal density matrices 
in the energy representation satisfy it. Therefore we cannot fix the Oth order term p from the Oth 
order equation (l29l ). 



To handle this problem, we introduce a projection superoperator P, which is defined by 

p\p\/P I _ / \ E i)( E i\ ( i= j) fm 

P\Eim\ - | {i¥:j) • (33) 

Namely, P is the projection to the diagonal part. We also define Q = 1 — P, which is the projection 
to the off-diagonal part. 

Hereafter we assume that H s is non-degenerate. Then the Oth order equation (1291 means that p is 
diagonal. Since C satisfies 

PC = C P = 0, (34) 

we obtain 

PC 1Po = PCiPpo = (35) 

from the 1st order equation d30l . Eq. (|35l ) means PC\P has a zero eigenvalue, and we assume that it 
is non-degenerate. Then Eq. (l35l > determines the Oth order term p uniquely. 

The unperturbed Liouvillian £ acts on the density matrix as 

(CoP), k = i(l»s,?]) jk (36) 

= ^j^P it , (37) 

in 

where Aj k = (Ei\A\E k ) denotes a matrix element in the energy representation. £ does not have 
its inverse because it has zero eigenvalues. Nevertheless we can define its inverse in the off-diagonal 
subspace: 

((QC Qr l p) jk = ^^Pik. (38) 

Then from (1301 we obtain the off-diagonal part of the first order term 

Q Pl = -(QCoQ)- 1 ^. (39) 



The diagonal part of the second order equation (131b can be rewritten as 

PdiP + Q)^ = 0. (40) 

P C\P has a zero eigenvalue, and the corresponding eigenvector is po- Therefore the general solution 
of ©is 

P Pl = -(P£ 1 P)'- 1 P£ 1 Q Pl + 1Po , (41) 

where (P£\P) _1 is the inverse of PC\P in the subspace spanned by non-zero eigenvectors of PC\P, 
and 7 is a number determined by the normalization condition Trpi = 0. 

In the same procedure we can determine Qp2, Pf>2, Qp3 and so on. These higher order terms, however, 
may be physically irrelevant because Eq. (|25l was derived from the approximation up to the first order 
off. 



3.2 Perturbative expansion with respect to A/3 

The pertubative solution we have obtained in the previous subsection is still very formal because we 
do not know the explicit form of the Oth order term p . In this and following subsections we try to 
find a more explicit solution by expanding p with respect to A/3, the inverse temperature difference 
between the two heat bath. 



We put 



0L = £-4p (42) 



Then we expand the heat bath superoperators and the density matrix as 



T L (/3 L ) = r L (/3)-^r L (/3) + 0(A/3 2 ), (44) 

T R (P R ) = T R ^) + ^T L ((3) + 0(A(3 2 ), (45) 

(46) 

p = p 00 + Af3p 01 +v(p 10 + Af3p u ) + O(v 2 ) + O(Af3 2 ). (47) 

We obtain an equation for each order 0(v n Af3 m ): 



0(1) 

0{A(5) 

Oh) 



£oPoo = 0, (48) 

£oPoi = 0, (49) 

C 0Pl0 + (T L + T R )p 00 = 0, (50) 



0(vA(3) : Copu + (T L + T R )p m + ^{-d p T L + d p Y R )p m = 0. (51) 

Note that T a and dpT a in the above equations are evaluated at the inverse temperature /3. 
Eqs. (l48l) and (l49l) mean that p 00 and p i are diagonal. Since T% satisfies 

PT%P = 0, (52) 

we obtain 

P(it + rf)/3oo = (53) 

by applying P to (|5Q|) . It has the equilibrium solution 

Poo = |e"^ s , (54) 



where Z is the partition function. Then from (|501 we obtain 



Qp 10 = -^{QC,Q)-\T L 2 +T R )e-^ s . (55) 



TlH s Tl = 


~- H s , 


UXj'U = 




YiRffi - 


hR 


nwjjn = 


= &$ 



3.3 Symmetric case 

By applying P to (I5TT) we obtain 

P(T L X + rf )p i + \p(-dpT{ + fyrf )p o = 0. (56) 

In principle, p i is determined by solving this equation. However, the inverse of -P(Tf + rf )P is hard 
to calculate analytically. 

Here we assume that the system and the heat baths are reflection symmetric. More precisely, we 
assume that YiH toX t[ = -f^tot, where ti is the parity operator which satisfies ft 2 = 1. Then we have 

(57) 
(58) 
(59) 
(60) 

Note that the operators are evaluated at the same inverse temperature (3 in Eqs. (1591 and (1601) . 
Then let us consider the second term of Eq. (l56t . Since p 00 is symmetric, we have 

n^rfpoon = ^rfp 00 (6i) 

n^rfpoofi = fyrfpoo. (62) 

A diagonal element in the second term of (l56l) is 

(E p \(-d^ + d^)p 00 \E p ) (E p \tl(-d^ + d Tf)p oo tl\E p ) (63) 

tl\E p ) = ±|£ p }) (64) 

(E p \(-d fi rf + d p rf)poo\E p ) (65) 

= -(E p \(-d^ + d^)p 00 \E p ). (66) 

Hence 

(E p \{-d p rf + d p rf)poo\E p ) = (67) 

and the second term of (l56l) vanishes. Then we have 

p(rf + rf )A)i = o, (68) 

whose solution is 

Poi oc p 00 = i e -' 3 ^. (69) 

To keep the normalization condition Trp = 1, we should put 

Poi = 0. (70) 

Then we obtain the lowest order nonequilibrium term 

Qpn = -^(QC,Q)- l Q(-dpT L + dpT R )e-^ s . (71) 

from (|5TT >. This is our main result. Note that diagonal elements do not contribute to nonequilibrium 
properties like the energy current and the temperature gradient. 



4 Energy current 

4.1 Energy current operator 

Let us consider the energy current going through the system. We divide the system into two parts. 
Then the system Hamiltonian is 

H s = H t + H % + H r , (72) 

where Hi and H r are the Hamiltonians for the left and right parts of the system, respectively, and Hi 
is the interaction between them. Note that [Hi, H r ] = 0. The energy current which goes from the left 
part to the right part can be defined as the energy loss of the left part: 

Ji = -H l = ~[Hi,H S ] = ~[Hi,H i ]. (73) 

in in 

It is also possible to define another current operator by the energy gain of the right part: 

I = H r = hH r ,H s ] = ^-[H r ,H t ]. (74) 

in in 

We can also define the energy current at a boundary between the system and a heat bath. The total 
energy of the system changes as 

±<Hs) = Tr(^p) (75) 

= Tr{H s ([H 3 ,p]+vr L p + vr R p}} (76) 

= wTr (# s r L p) + vTr (^H s r R pj . (77) 

The left (right) term can be interpreted as the energy current at the left (right) boundary. Therefore 
we define two current operators Jl and Jr so that the following relations hold. 

(Jl) = vTr (H s r L p), (78) 

(J R ) = -vTr (h s T r p) . (79) 

Then 



Tr(J L p) (80) 

-^2 E Tr fa {[xf, RfiP] + [xj, Rf lP ] ] ) + i {[Xf, Wjip] - [Xf, Wffi) } (81) 
jl 

J^Tr | ([H s ,Xf] (AJ + iWJi) + [k% + iW$f [H S ,X^ p\ . (82) 



2h 2 
jl 



Hence 



Jl = -^E{[^^K^ + ^) + (^ + ^) t[ ^ 5 '^ ]t }- (83) 



jl 

In the same way we obtain 



J * = ^E{[^^f](4l + ^) + (^ + ^) t fe^f] t }- (84) 

ji 

In the steady state all current operators should have the same expectation value. 



4.2 Expectation value 

Let us consider the expectation value of a energy current operator in the system. In our perturbation 
theory, Qpn is the leading nonequilibrium term. Therefore we evaluate 

Tr(JiQpu) = ^{[Hi,H s ](QC Q)-\-d^ L + d^ R )p 00 } (85) 

= —[{CoH^QCoQ^Qi-d^ + dpF^pooY (86) 

Note that 

Tr {(£ A) (QC^Qy'QB} = -Tr(AQB) (87) 

holds in general, because 

Tr {(C A)(QC; 1 Q)- 1 QB'j = ^E^C^E^^E^QC^Q^QB^) (88) 

l,m 

= y j (E l -E m )(E l \A\E m }—^—-(E m \B\E l } (89) 

' Elm. — -til 






-J]<£,|i|£ ro )<£ ro |£|£,) (90) 

-Tr(AQB). (91) 



Hence 



Tr^Qpn) = ^Tr {#,Q(-fyr L + d r R )p m } (92) 

= iTr{^(-^r L H-^r R )p 00 } (93) 

(■.■P(-dpr L + dpY R )p 00 = 0). (94) 

Each term in dpT a p m has the form [1", dpZjipw] (Z = R, W). Then 

T* (£,[*?, fy^poo]) = Tr([^,X;]^/)oo), (95) 
which vanishes if a = /?. Hence 

TriJiQpu) = -\Tr (H l d fi r L p 0Q ^ (96) 

= -i-Tr (ffsfyl^poo) (97) 

= -^TV [h s (d $ V{ + d F%) poo}- (98) 
In the second line we used [Hi, Xj] = [Hs, Xj}. With some algebra, one can easily show 

Tr (ifsfyltpoo) = 0. (99) 
Since 

Tf (/3)e-^ s = (100) 



for any f3, 

Op (t l x (/3)e-^ s ) = dpT[ (/3)e-^ s + T[ (/3)^e"^ = 0. (101) 

Therefore 

^rfpoo = -T^HsPoo (102) 



and 



TrUQpn) = \tt (h S ^HsPoo) (103) 

l -Re J2 Tr (H s [Xf, R^H sPm ]) (104) 



2ft 2 



^Re J2 Tr ([#s, Xfl^^poo) (105) 



2ft 2 



^ E E {^(^ - E q )(E p \X}\E q )(E q \Xt\E p )% 1 (u qp )e-^ E ' 



Ah 2 Z 

ji VI 



+ E P (E P - E q )(E q \Xf\E p )(E p \X^\E q )%(u qp )e-^ (106) 

^ E E { E p( E p - E ,)( E p\X}\ E <H E <\X?\ E p)*riuv)e- p * 



4h 2 Z 

ji pq 



+ E q (E q - E p ){E p \Xf\E q ){E q \X^\E p )^{uj pq )e-^ (107) 

= -4^EE{^- B «) ! ^M^WW^} (io8) 

= -i^E(^'^H^,^]) /3 . (109) 

J' 

In the last line, the expectation value is evaluated for the system in equilibrium at the inverse temper- 
ature (5. 

In the same way we obtain 

Tr(J r g Pll ) = --^Y,([Hs,X*\[R*,H s ]) . (HO) 

ji 

Since we have assumed the reflection symmetry, we obtain 

(Ji) = (l) = -^zZ(^Xf][R^,H s }) vA(3 + 0(v 2 ) + 0(A(3 2 ). (Ill) 



We can also calculate the expectation values of the current operators at the boundaries. Let us consider 
Jl- Since Jl is 0(v), the lowest order current comes from O(v ) terms of the density matrix. Because 
Poi = in the symmetric case, we have 

(J L ) = Tr(J L p 00 ) + O(v 2 ) + O(A/3 2 ). (112) 



Then, from Eq. (|78T ) 

Tr(J L p 00 ) = vTr(H s r L p oo y (113) 

Note that T L is evaluated at (3 L = (3 - A/3/2 here. Substituting 

TlWl) = T L {p) - ^-d p T L (P) (114) 

we obtain 

Tr(J L p 00 ) = "^Tr (H s d p Y L p m ) , (115) 

which is equivalent to (|97l ). Then we obtain the same current expectation value again: 

(Jl) = (Jr) = -^J2(fts,X?][R?i,Hs}) vA(3 + 0(v 2 ) + 0(Af3 2 ). (116) 

ji 



4.3 Kubo formula 

Let us consider the total system including the heat bath again. The current operator for the energy 
coming from the heat bath L to the system is defined as 

J'l = j t (H s + uH RS + H R ) (117) 

= ±_[H s + uH RS + H R ,H tot } (118) 

= ^[H s ,uH LS ] (119) 



ih 

3 



Then we define the correlation function 



C L (t) = \{J' L {t)J'L + J'J'L{t))- (121) 

The expectation value is evaluated for the equilibrium of the total system with inverse temperature (3. 
Since J' L is 0(u 2 ), we have 

C L (t) = C^(t) + 0(u 3 ), (122) 

where 

C (0) (t) = -Tr {pf <g) pg f e KHs+H B )t/nj> Le -i(H s+ H B )t/nj, L + j^As+AbWkj^&s+AbW^ j (123) 

Here, p^ q and jf^ represents the equilibrium state of the system and the heat baths, respectively. Then 
we have 






►124) 



Jo 



and 

i>oo 2 

dtCf\t) = ~EE^(^-£ f )W|^)(^) 

/■OO 

x / di{e _iw «'*$ i ,(*) + e*' , « ,t $ i i(-i)} (125) 

= -^EE e ^ Ep ^-^) 2 (^l^il^}(^l^l^}^K P ) (126) 



v 



J2([Hs,Xf}[Rf,H s }\ (127) 



2ft 2 
ji 

Therefore, in the lowest order, the current is written as 

(Jl) = ^ J dtC L (t), (128) 

which is the Kubo formula [3], or the fluctuation-dissipation theorem, in this case. Note that A/3/2 is 
the inverse temperature difference at the boundary. 

In spite of the formal similarity, physical content of Eq. (11281) is quite different from the original 
Kubo formula [1]. For example, the transport coefficient contains the information of the heat baths, 
though the original one does not. 



5 Summary 

We have calculated the density matrix for the NESS using the time-independent perturbation theory. 
Our main result is Eq. (1711 . which is an explicit expression for the density matrix for the NESS in the 
reflection symmetric setting. We have also calculated the expectation value of the energy current and 
shown that the Kubo formula holds in this case. 
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